#include "DXUT.h"
#include "PreethamSkyLight.h"

/*
	File:		SunSky.cc
 
    Function:	Implements a sky radiance function for a clear sky
 
    Author:     Brian Smits (bes@phoenix.cs.utah.edu),
    			Some minor modifications by Andrew Willmott (ajw@cs.cmu.edu)
 
	Notes:		
*/

/**********************************************************
//  TODO							   
// Fix Indexing (make it zero based, not 1 based...)
// Remove degrees and radians macros
// Define SUN_CONSTANT
// Confirm I am calculating Y correctly for RiColorXYZV
// ********************************************************/


/**********************************************************
// Atmospheric Perspective Constants
// 
// ********************************************************/

// Aerial Perspective
// Oct 23, 1998: Started
// Nov  3, 1998: Switched from 1-particle to 2-particle model.
// Nov 13, 1998: Inheritance.
// Dec  4, 1998: Support for turbidity > 2
// Dec  9, 1998: Code under ggSunSky
// Dec 11, 1998: Support for floating turbidity
// Feb 22, 1999: 3-particle (fog) model.

// density decreses as exp(-alpha*h) alpha in m^-1.
// _1 refers to Haze, _2 referes to molecules.


static float const Alpha_1 = 0.83331810215394e-3;
static float const Alpha_2 = 0.11360016092149e-3;

const int kThetaBins = 20;
const int kPhiBins = 20;

static Spectrum AT0_1[kThetaBins + 1][kPhiBins + 1];
static Spectrum AT0_2[kThetaBins + 1][kPhiBins + 1];

inline Spectrum exp(const Spectrum &in)
{
	Spectrum out;
	for (Int i = 0; i < Spectrum::kComponents; i++)
		out[i] = exp(in[i]);
	return out;
}

/* All angles in radians, theta angles measured from normal */
inline float AngleBetween(float thetav, float phiv, float theta, float phi)
{
	float cospsi = sin(thetav) * sin(theta) * cos(phi - phiv) + cos(thetav) * cos(theta);

	if (cospsi > 1.0)
		return 0;
	if (cospsi < -1.0)
		return vl_pi;

	return acos(cospsi);
}

PreethamSkyLight::PreethamSkyLight(
		    GCLReal		lat,
		    GCLReal 	longi,
		    int 		sm,        		// standardMeridian
		    int 		jd,        		// julianDay
		    GCLReal		tOfDay,    		// timeOfDay
		    GCLReal 	turb, 			// turbidity
		    bool		initAtmEffects
)
{
	latitude = lat;
	longitude = longi;
	julianDay = jd;
	timeOfDay = tOfDay;
	standardMeridian = int(sm * 15.0);   // sm is actually timezone number (east to west, zero based...)
	turbidity = turb;

	V = 4.0;  // Junge's exponent.

	InitSunThetaPhi();
	toSun = D3DXVECTOR3(cos(phiS) * sin(thetaS), sin(phiS) * sin(thetaS), cos(thetaS));

	// Units: W.cm^-2.Sr^-1.
	sunSpectralRad = ComputeAttenuatedSunlight(thetaS, turbidity);
	sunSolidAngle = 0.25 * vl_pi * 1.39 * 1.39 / (150 * 150);   // = 6.7443e-05

	GCLReal theta2 = sqr(thetaS);
	GCLReal theta3 = theta2 * thetaS;
	GCLReal T = turb;
	GCLReal T2 = sqr(turb);

	GCLReal chi = (4.0 / 9.0 - T / 120.0) * (vl_pi - 2 * thetaS);
	zenith_Y = (4.0453 * T - 4.9710) * tan(chi) - .2155 * T + 2.4192;
	zenith_Y *= 1000;   // conversion from kcd/m^2 to cd/m^2

	zenith_x =
	    ( + 0.00165 * theta3 - 0.00374 * theta2 + 0.00208 * thetaS + 0) * T2 +
	    ( -0.02902 * theta3 + 0.06377 * theta2 - 0.03202 * thetaS + 0.00394) * T +
	    ( + 0.11693 * theta3 - 0.21196 * theta2 + 0.06052 * thetaS + 0.25885);

	zenith_y =
	    ( + 0.00275 * theta3 - 0.00610 * theta2 + 0.00316 * thetaS + 0) * T2 +
	    ( -0.04214 * theta3 + 0.08970 * theta2 - 0.04153 * thetaS + 0.00515) * T +
	    ( + 0.15346 * theta3 - 0.26756 * theta2 + 0.06669 * thetaS + 0.26688);

	perez_Y[0] = 0.17872 * T - 1.46303;
	perez_Y[1] = -0.35540 * T + 0.42749;
	perez_Y[2] = -0.02266 * T + 5.32505;
	perez_Y[3] = 0.12064 * T - 2.57705;
	perez_Y[4] = -0.06696 * T + 0.37027;

	perez_x[0] = -0.01925 * T - 0.25922;
	perez_x[1] = -0.06651 * T + 0.00081;
	perez_x[2] = -0.00041 * T + 0.21247;
	perez_x[3] = -0.06409 * T - 0.89887;
	perez_x[4] = -0.00325 * T + 0.04517;

	perez_y[0] = -0.01669 * T - 0.26078;
	perez_y[1] = -0.09495 * T + 0.00921;
	perez_y[2] = -0.00792 * T + 0.21023;
	perez_y[3] = -0.04405 * T - 1.65369;
	perez_y[4] = -0.01092 * T + 0.05291;

	if (initAtmEffects)
	{
		// Atmospheric Perspective initialization
		CreateConstants();
		InitA0();
	}

	atmInited = initAtmEffects;
}

/**********************************************************
// South = x,  East = y, up = z
// All times in decimal form (6.25 = 6:15 AM)
// All angles in Radians
// ********************************************************/

void PreethamSkyLight::GetSunThetaPhi(GCLReal &stheta, GCLReal &sphi) const
{
	sphi = phiS;
	stheta = thetaS;
}

/**********************************************************
// South = x,  East = y, up = z
// All times in decimal form (6.25 = 6:15 AM)
// All angles in Radians
// ********************************************************/
void PreethamSkyLight::InitSunThetaPhi()
{
	GCLReal solarTime = timeOfDay +
	                    (0.170 * sin(4 * vl_pi * (julianDay - 80) / 373) - 0.129 * sin(2 * vl_pi * (julianDay - 8) / 355)) +
	                    (standardMeridian - longitude) / 15.0;

	GCLReal solarDeclination = (0.4093 * sin(2 * vl_pi * (julianDay - 81) / 368));

	GCLReal solarAltitude = asin(sin(DegsToRads(latitude)) * sin(solarDeclination) -
	                             cos(DegsToRads(latitude)) * cos(solarDeclination) * cos(vl_pi * solarTime / 12));

	GCLReal opp, adj;
	opp = -cos(solarDeclination) * sin(vl_pi * solarTime / 12);
	adj = -(cos(DegsToRads(latitude)) * sin(solarDeclination) +
	        sin(DegsToRads(latitude)) * cos(solarDeclination) * cos(vl_pi * solarTime / 12));
	GCLReal solarAzimuth = atan2(opp, adj);

	phiS = -solarAzimuth;
	thetaS = vl_pi / 2.0 - solarAltitude;
}

Vector PreethamSkyLight::GetSunPosition() const
{
	return toSun;
}

/**********************************************************
// Solar Radiance
// 
// ********************************************************/
Spectrum PreethamSkyLight::GetSunSpectralRadiance() const
{
	return sunSpectralRad;
}

GCLReal PreethamSkyLight::GetSunSolidAngle() const
{
	return sunSolidAngle;
}

/**********************************************************
// Sky Radiance
// 
// ********************************************************/
Spectrum PreethamSkyLight::GetSkySpectralRadiance(const Vector &varg) const
{
	GCLReal	theta, phi;
	Vector	v = norm(varg);

	if (v[2] < 0) return Spectrum(0.000000);
	if (v[2] < 0.001)
		v = norm(Vector(v[0], v[1], 0.001));

	theta = acos(v[2]);
	if (fabs(theta) < 1e-5)
		phi = 0;
	else
		phi = atan2(v[1], v[0]);

	return GetSkySpectralRadiance(theta, phi);
}

inline GCLReal PreethamSkyLight::PerezFunction(const GCLReal *lam, GCLReal theta, GCLReal gamma, GCLReal lvz) const
{
	GCLReal den = ((1 + lam[0] * exp(lam[1])) *
	               (1 + lam[2] * exp(lam[3] * thetaS) + lam[4] * cos(thetaS) * cos(thetaS)));

	GCLReal num = ((1 + lam[0] * exp(lam[1] / cos(theta))) *
	               (1 + lam[2] * exp(lam[3] * gamma) + lam[4] * cos(gamma) * cos(gamma)));

	return(lvz * num / den);
}

Spectrum PreethamSkyLight::GetSkySpectralRadiance(GCLReal theta, GCLReal phi) const
{
	GCLReal gamma = AngleBetween(theta, phi, thetaS, phiS);
	// Compute xyY values
	GCLReal x = PerezFunction(perez_x, theta, gamma, zenith_x);
	GCLReal y = PerezFunction(perez_y, theta, gamma, zenith_y);
	GCLReal Y = PerezFunction(perez_Y, theta, gamma, zenith_Y);

	Spectrum spect = ChromaticitySpectrum(x, y);

	return((Y * spect) / spect.ToCIE_XYZ()[1]);
}

Colour PreethamSkyLight::GetSkyXYZRadiance(const Vector &varg) const
{
	GCLReal	theta, phi;
	Vector	v = norm(varg);

	if (v[2] < 0) return cBlack;
	if (v[2] < 0.001)
		v = norm(Vector(v[0], v[1], 0.001));

	theta = acos(v[2]);
	if (fabs(theta) < 1e-5)
		phi = 0;
	else
		phi = atan2(v[1], v[0]);

	return(GetSkyXYZRadiance(theta, phi));
}

Colour PreethamSkyLight::GetSkyXYZRadiance(GCLReal theta, GCLReal phi) const
{
	GCLReal gamma = AngleBetween(theta, phi, thetaS, phiS);
	// Compute xyY values
	GCLReal x = PerezFunction(perez_x, theta, gamma, zenith_x);
	GCLReal y = PerezFunction(perez_y, theta, gamma, zenith_y);
	GCLReal Y = PerezFunction(perez_Y, theta, gamma, zenith_Y);

	GCLReal X = (x / y) * Y;
	GCLReal Z = ((1 - x - y) / y) * Y;

	return(Colour(X, Y, Z));
}


/**********************************************************
// Atmospheric perspective functions
// 
// ********************************************************/


/**********************************************************
// Initialization
// 
// ********************************************************/
void PreethamSkyLight::CalculateA0(GCLReal thetav, GCLReal phiv, Spectrum &A0_1, Spectrum &A0_2) const
{
	GCLReal psi;
	Spectrum skyRad;
	Spectrum beta_ang_1, beta_ang_2;
	GCLReal theta, phi, phiDelta = vl_pi / 20;
	GCLReal thetaDelta = vl_pi / 2 / 20;
	GCLReal thetaUpper;

	Spectrum skyAmb_1 = 0;
	Spectrum skyAmb_2 = 0;

	thetaUpper = vl_pi / 2;

	for (theta = 0; theta < thetaUpper; theta += thetaDelta)
		for (phi = 0; phi < 2 * vl_pi; phi += phiDelta)
		{
			skyRad = GetSkySpectralRadiance(theta, phi);
			psi = AngleBetween(thetav, phiv, theta, phi);

			beta_ang_1 = beta_p_ang_prefix * GetNeta(psi, V);
			beta_ang_2 = beta_m_ang_prefix * (1 + 0.9324 * cos(psi) * cos(psi));

			skyAmb_1 += skyRad * beta_ang_1 * sin(theta) * thetaDelta * phiDelta;
			skyAmb_2 += skyRad * beta_ang_2 * sin(theta) * thetaDelta * phiDelta;
		}
	/* Sun's ambience term*/
	psi = AngleBetween(thetav, phiv, thetaS, phiS);

	beta_ang_1 = beta_p_ang_prefix * GetNeta(psi, V);
	beta_ang_2 = beta_m_ang_prefix * (1 + 0.9324 * cos(psi) * cos(psi));

	Spectrum sunAmb_1 = sunSpectralRad * beta_ang_1 * sunSolidAngle;
	Spectrum sunAmb_2 = sunSpectralRad * beta_ang_2 * sunSolidAngle;

	// Sum of sun and sky  (should probably add a ground ambient)
	A0_1 = (sunAmb_1 + skyAmb_1);
	A0_2 = (sunAmb_2 + skyAmb_2);
}

Void PreethamSkyLight::InitA0() const
{
	int i, j;
	GCLReal theta, phi;

	GCLReal delTheta = vl_pi / kThetaBins;
	GCLReal delPhi = 2.0 * vl_pi / kPhiBins;

	cerr << "SunSky::Preprocessing: 0%\r";
	for (i = 0, theta = 0; theta <= vl_pi; theta += delTheta, i++)
	{
		for (j = 0, phi = 0; phi <= 2.0 * vl_pi; phi += delPhi, j++)
			CalculateA0(theta, phi, AT0_1[i][j], AT0_2[i][j]);
		cerr << "SunSky::Preprocessing: " << 100 * (i * kPhiBins + j)
		/ ((kThetaBins + 1) * kPhiBins) << "%  \r";
	}
	cerr << "SunSky::Preprocessing:  100%   " << endl;
}

/**********************************************************
// Evaluation
// 
// ********************************************************/
void PreethamSkyLight::GetAtmosphericEffects(const Vector &viewer, const Vector &source,
                                   Spectrum &attenuation, Spectrum &inscatter ) const
{
	Assert(atmInited, "atm not initted");
	//#define LARGE_DISTANCE 300000
	// Clean up the 1000 problem
	GCLReal h0 = viewer[2] + 1000;   //1000 added to make sure ray doesn't
	//go below zero.
	Vector direction = norm(source - viewer);
	GCLReal thetav = acos(direction[2]);
	GCLReal phiv = atan2(direction[1], direction[0]);
	GCLReal s = len(viewer - source);

	// This should be changed so that we don't need to worry about it.
	if (h0 + s * cos(thetav) <= 0)
	{
		attenuation = 1;
		inscatter = 0;
		cerr << "\nWarning: Ray going below earth's surface \n";
		return;
	}

	attenuation = AttenuationFactor(h0, thetav, s);
	inscatter = InscatteredRadiance(h0, thetav, phiv, s);
}


inline GCLReal EvalFunc(GCLReal B, GCLReal x)
{
	if (fabs(B*x) < 0.01) return x;
	return (1 - exp( -B*x)) / B;
}


Spectrum PreethamSkyLight::AttenuationFactor(GCLReal h0, GCLReal theta, GCLReal s) const
{
	GCLReal costheta = cos(theta);
	GCLReal B_1 = Alpha_1 * costheta;
	GCLReal B_2 = Alpha_2 * costheta;
	GCLReal constTerm_1 = exp( -Alpha_1 * h0) * EvalFunc(B_1, s);
	GCLReal constTerm_2 = exp( -Alpha_2 * h0) * EvalFunc(B_2, s);

	return (exp( -beta_p * constTerm_1) *
	        exp( -beta_m * constTerm_2));
}


static void GetA0fromTable(GCLReal theta, GCLReal phi, Spectrum &A0_1, Spectrum &A0_2)
{
	GCLReal eps = 1e-4;

	if (phi < 0) phi += 2 * vl_pi;  // convert phi from -pi..pi to 0..2pi
	theta = theta * kThetaBins / vl_pi - eps;
	phi = phi * kPhiBins / (2 * vl_pi) - eps;
	if (theta < 0) theta = 0;
	if (phi < 0) phi = 0;
	Int i = (Int) theta;
	GCLReal u = theta - i;
	Int j = (Int)phi;
	GCLReal v = phi - j;

	A0_1 = (1 - u) * (1 - v) * AT0_1[i][j] + u * (1 - v) * AT0_1[i + 1][j]
	       + (1 - u) * v * AT0_1[i][j + 1] + u * v * AT0_1[i + 1][j + 1];
	A0_2 = (1 - u) * (1 - v) * AT0_2[i][j] + u * (1 - v) * AT0_2[i + 1][j]
	       + (1 - u) * v * AT0_2[i][j + 1] + u * v * AT0_2[i + 1][j + 1];
}

inline GCLReal Helper1(GCLReal A, GCLReal B, GCLReal C, GCLReal D,
                         GCLReal H, GCLReal K, GCLReal u)
{
	GCLReal t = exp( -K * (H - u));
	return (t / K)*((A*u*u*u + B*u*u + C*u + D) -
	                (3*A*u*u + 2*B*u + C) / K +
	                (6*A*u + 2*B) / (K*K) -
	                (6*A) / (K*K*K));
}

inline void CalculateABCD(GCLReal a, GCLReal b, GCLReal c, GCLReal d, GCLReal e,
                            GCLReal den, GCLReal &A, GCLReal &B, GCLReal &C, GCLReal &D)
{
	A = ( -b * d - 2 + 2 * c + a * e - b * e + a * d) / den;
	B = -(2 * a * a * e + a * a * d - 3 * a - a * b * e +
	      3 * a * c + a * b * d - 2 * b * b * d - 3 * b - b * b * e + 3 * b * c) / den;
	C = ( -b * b * b * d - 2 * b * b * a * e - b * b * a * d + a * a * b * e +
	      2 * a * a * b * d - 6 * a * b + 6 * b * a * c + a * a * a * e) / den;
	D = -( b * b * b - b * b * b * a * d - b * b * a * a * e + b * b * a * a * d
	       - 3 * a * b * b + b * e * a * a * a - c * a * a * a + 3 * c * b * a * a) / den;
}

Spectrum PreethamSkyLight::InscatteredRadiance(GCLReal h0, GCLReal theta, GCLReal phi, GCLReal s) const
{
	GCLReal costheta = cos(theta);
	GCLReal B_1 = Alpha_1 * costheta;
	GCLReal B_2 = Alpha_2 * costheta;
	Spectrum A0_1;
	Spectrum A0_2;
	Spectrum result;
	Int i;
	Spectrum I_1, I_2;

	GetA0fromTable(theta, phi, A0_1, A0_2);

	// approximation (< 0.3 for 1% accuracy)
	if (fabs(B_1*s) < 0.3)
	{
		GCLReal constTerm1 = exp( -Alpha_1 * h0);
		GCLReal constTerm2 = exp( -Alpha_2 * h0);

		Spectrum C_1 = beta_p * constTerm1;
		Spectrum C_2 = beta_m * constTerm2;

		for (Int i = 0; i < Spectrum::kComponents; i++)
		{
			I_1[i] = (1 - exp( -(B_1 + C_1[i] + C_2[i]) * s)) / (B_1 + C_1[i] + C_2[i]);
			I_2[i] = (1 - exp( -(B_2 + C_1[i] + C_2[i]) * s)) / (B_2 + C_1[i] + C_2[i]);
		}

		return A0_1 * constTerm1 * I_1 + A0_2 * constTerm2 * I_2;
	}

	// Analytical approximation
	GCLReal A, B, C, D, H1, H2, K;
	GCLReal u_f1, u_i1, u_f2, u_i2, Int_f, Int_i, fs, fdashs, fdash0;
	GCLReal a1, b1, a2, b2;
	GCLReal den1, den2;

	b1 = u_f1 = exp( -Alpha_1 * (h0 + s * costheta));
	H1 = a1 = u_i1 = exp( -Alpha_1 * h0);
	b2 = u_f2 = exp( -Alpha_2 * (h0 + s * costheta));
	H2 = a2 = u_i2 = exp( -Alpha_2 * h0);
	den1 = (a1 - b1) * (a1 - b1) * (a1 - b1);
	den2 = (a2 - b2) * (a2 - b2) * (a2 - b2);

	for (i = 0; i < Spectrum::kComponents; i++)
	{
		// for Integral 1
		K = beta_p[i] / (Alpha_1 * costheta);
		fdash0 = -beta_m[i] * H2;
		fs = exp( -beta_m[i] / (Alpha_2 * costheta) * (u_i2 - u_f2));
		fdashs = -fs * beta_m[i] * u_f2;

		CalculateABCD(a1, b1, fs, fdash0, fdashs, den1, A, B, C, D);
		Int_f = Helper1(A, B, C, D, H1, K, u_f1);
		Int_i = Helper1(A, B, C, D, H1, K, u_i1);
		I_1[i] = (Int_f - Int_i) / ( -Alpha_1 * costheta);

		// for Integral 2
		K = beta_m[i] / (Alpha_2 * costheta);
		fdash0 = -beta_p[i] * H1;
		fs = exp( -beta_p[i] / (Alpha_1 * costheta) * (u_i1 - u_f1));
		fdashs = -fs * beta_p[i] * u_f1;

		CalculateABCD(a2, b2, fs, fdash0, fdashs, den2, A, B, C, D);
		Int_f = Helper1(A, B, C, D, H2, K, u_f2);
		Int_i = Helper1(A, B, C, D, H2, K, u_i2);
		I_2[i] = (Int_f - Int_i) / ( -Alpha_2 * costheta);

	}
	return A0_1 * I_1 + A0_2 * I_2;
}
